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Abstract 

Studies of the many-body potential surface of liquid sodium have 



■ shown that it consists of a great many intersecting nearly harmonic 

, valleys, a large fraction of which have the same frequency spectra. 

| This suggests that a sufficiently supercooled state of this system, re- 

• maining in a single valley, would execute nearly harmonic motion. To 

test this hypothesis, we have compared Z(t), the normalized veloc- 
ity autocorrelation function, calculated from MD simulations to that 
predicted under the assumption of purely harmonic motion. We find 
nearly perfect agreement between the two, suggesting that the har- 
monic approximation captures all essential features of the motion. 

G 

O 

P. : 1 Introduction 

> ■ 

^ ■ Recent work by Wallace and Clements [I], [| has uncovered several important 

properties of the many-body potential underlying the motion of liquid sodium 
systems. Specifically, it has been shown that (a) the potential surface consists 
of a large number of intersecting nearly harmonic valleys, (b) these valleys 
can be classified as symmetric (crystalline, micro crystalline, or retaining some 
nearest-neighbor remnants of crystal symmetry) or random, with the random 
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valleys vastly outnumbering the symmetric ones, (c) the frequency spectra 
of different random valleys are nearly identical (while those of the symmetric 
valleys vary widely), and (d) below 35 K the system remains in a single valley 
throughout the longest molecular dynamics (MD) runs that were performed. 
Results (a) through (c) verify predictions made by Wallace in his theory of 
liquid dynamics ||, which has been successfully applied to account for the 
high-temperature specific heats of monatomic liquids [|J and a study of the 
velocity autocorrelation function || . These four results together suggest that 
below 35 K the motion of the atoms in liquid sodium is purely harmonic to 
a high degree of approximation, again as predicted by Wallace in [[J , and we 
would like to test this hypothesis further. One check is to compare the mean 
square displacement from MD with the prediction from purely harmonic 
motion, which is done in Fig. 12 of [J], where the two are found to agree 
closely. However, it would be more convincing if the theory could be shown 
to reproduce an entire scalar function calculated from MD (instead of just 
a single number), such as the normalized velocity autocorrelation function 
Z(t). That is the aim of this paper. We will show that purely harmonic 
motion of the atoms in a potential valley produces a Z(t) which matches 
that of MD calculations to within the calculations' accuracy; thus we will 
conclude that the motion of atoms in a nondiffusing supercooled liquid state 
is very nearly entirely harmonic. For completeness, in Sec. |2| we briefly review 
the calculation of Z(t) assuming harmonic motion, and in Sec. |3| we compare 
this result with MD. Finally, in Sec. [|we make contact with work by others 
in this field, as well as comparing these results to Wallace's earlier effort || 
mentioned above. 



2 Harmonic Theory 

If an iV-body system is moving in a potential valley, the potential can be 
expanded about the valley minimum with the resulting Hamiltonian 

H = £' fjjT + E ' ^LjUKiULj + $ A (1) 
Ki Z1V1 Ki,Lj 

where uk% is the ith component of the Kth particle's displacement from 
equilibrium, pxi is the corresponding momentum, and the anharmonic term 
&a contains all of the higher order parts of the expansion. The primed sum 
indicates that the sum is performed under the constraint that the center of 
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mass of the system is stationary. (As a result, the system has only 3iV — 3 
independent degrees of freedom.) The matrix &Ki,Lj is called the dynamical 
matrix of the system. If the valley is approximately harmonic, we can neglect 
If coordinates q\ are defined by the relation 

u Ki = W Ki,\q\ (2) 
A 

where the Wxi,\ form a 3iV x 3iV orthogonal matrix, satisfying 

J2 W Ki,\ w Ki,\> = 5\\>, (3) 
Ki 

then the Hamiltonian in these new coordinates is 

H = J2wb+ 2 J2 W Ki^Ki,LjW Lj ,x>q\q\> (4) 
A Ki,Lj AA' 

where the p x are the momenta conjugate to the q x . Now one can always 
choose the wk%,x to diagonalize $Ki,Lj, so that 

w Ki,\^Ki,LjW L j,\' = Mu\8x\>. (5) 

Ki,Lj 

(This equation defines the frequencies u>\ in terms of the eigenvalues of 
®Ki,Lj-) With this choice, the Hamiltonian becomes 

" - ? (# + \ M ^) ■ < 6 » 

Three of the uj\ are zero; these modes correspond to uniform motion of the 
center of mass. Since we have restricted the center of mass position and 
velocity to zero, these modes are not excited. The classical equations of 
motion for the remaining modes are solved by 

q\(t) = a x sm(uj x t + «a), (7) 

or, returning to the original coordinates, 

UKi(t) = w Ki,\a\ sm(uj x t + a x ), (8) 

A 



3 



with the understanding that the sum on A ranges from 1 to 3iV — 3. The 
velocities of the particles are 

v Ki(t) = Y w Ki,x u x a x cos(uj\t + a A ). (9) 

A 

We compute the (v(t) ■ v(0)) in Z(t) by calculating v K (t) ■ v K (0), summing 
over K and dividing by N — 1 (remember that only 3iV — 3 coordinates are 
independent), and averaging over the amplitudes a x and phases a x . Thus 

Z(t) = \(v(t)-v(0)) 

Y Y w Ki,\w KitX >uJxuJx> (a x a x >)(cos(u x t + a x ) cos(a A ')) 



3^ - 3 Kl xx 



Y 8 xx >uj x ujy (a x a x ,)(cos(u x t + a x ) cos(a A ')) 



3^ - 3 AA , 



Y ^l^l) (cos(u x t + a x ) cos(a A )) 



SiV-3 A 



: Y"l( a l) cos M- ( 10 ) 



6N-6 A 



By the equipartition theorem, 



^) = \kT (11) 



for any nonzero cu x , from which it follows that 

2kT 



(a 2 x ) = (12) 
so 

1 kT 

Z(t) = oAr V cos(cu A t). (13) 

Notice that Z(0) = kT/M, so Z(t) defined by 

Z(t) = Z(0)Z(t) (14) 
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is given in this theory by 

^(0 = ^^5>os(u, A t). (15) 

A 

This is the result we wish to compare with MD. 

To do so, we need the frequencies ux, which are related to the eigenvalues 
of the dynamical matrix &Ki,Lj as indicated in Eq. (Q). These were evaluated 
for five separate random valleys in Q by quenching all the way down to a 
valley minimum and diagonalizing &Ki,Lj there; as pointed out in Sec. [I], 
these eigenvalues were found to be independent of the specific random valley 
chosen. All five sets of eigenvalues are shown in Fig. 7 of [|TJ, and we picked 
one set at random to use in performing the sum in Eq. (|i5|); the other sets 
produce identical graphs of Z(t). 

We can also use the set of eigenvalues to reconstruct the density of fre- 
quencies g(uj); the results are shown in Fig. [I]. Note that we do not actually 
integrate over this g(uj) to evaluate Z(t) below; we directly sum over the 
given set of frequencies as indicated in Eq. ( |15"D . The Figure is provided only 
to convey a sense of the shape of the frequency distribution. Also note that 
this g{uj) is determined from fully mechanical considerations; as a result, it 
is not temperature-dependent as are the frequency spectra used in Instanta- 
neous Normal Mode (INM) studies §, [7|, g |, [10|. We will expand on this 
point in the Conclusion. 



3 Comparison with MD 

The MD setup used to calculate Z(t) to compare with Eq. (|15"D is essentially 
that described in : N particles interact through a potential that is known 
to reproduce accurately a wide variety of experimental properties of metallic 
sodium (see discussion in |J for details). The two significant changes are 
that we used iV = 500 for all runs and that the MD timestep was reduced 
to 5t = 0.2t*, where t* = 7.00 x 10~ 15 s is the natural timescale defined 
in [1J. (The system's mean vibrational period r = 27r/c<j rms , where the rms 
frequency u rms is calculated in Jjj], is approximately 300 St.) We cooled the 
sodium sample to 22.3 K and 6.69 K, and then we ran each at equilibrium 
to collect velocities Vxit) to be used to calculate Z{t) by the formula 

Z(t) = 4E^tEMHO- v K (t'). (16) 
ojv K n+i t , =0 
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Figure 1: g(u) constructed from one of the five sets of frequencies in Fig. 7 
of [p. This same set of frequencies is used to calculate Z(t) from Eq. (|I5]) . 
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We then divided by Z(0) to obtain Z(t). The number n was chosen as large 
as possible without running beyond the data calculated in the MD run. We 
know that at these temperatures the sodium is nondiffusing for two reasons: 
Both temperatures are below the 35 K threshold and Z(t) from either 
MD run (shown below in Fig. ^) integrates to zero, yielding zero diffusion 
coefficient. 

The formula above may fail to produce reliable values of Z(t) for three 
reasons. First, the number of data points in the time average may be too 
small; if the MD simulation is run out to time t max , then for a given value of t 
in Eq. ([U]), the maximum possible value of the upper limit n is t max — t. Thus 
we require t max >> t; we have chosen t max = 50, 000 timesteps and we have 
calculated Z(t) only to t = 1000. To ensure that this value of t max is large 
enough, we also performed MD runs out to 200,000 timesteps and calculated 
Z(t) from them; the differences from the 50,000 timestep result were of order 
10~ 3 . Hence we are confident that 50,000 timesteps is enough if we calculate 
Z(t) to only 1000 timesteps. Second, it is possible that reducing the timestep 
(thus increasing the accuracy of the simulation) might improve the accuracy 
of Z(t). To test this, we performed another MD run with 5t reduced to 0.05t*, 
keeping the "real" time of the run the same; this also produced differences 
in Z(t) of order 10~ 3 . Thus we are sure that our timestep is small enough. 
Finally, there is the possibility of finite size effects. Since the MD system has 
periodic boundary conditions, an acoustic wave sent out from the system at 
t = could propagate across the simulation region and return to its point of 
origin in a finite time, producing spurious correlations that would show up 
in Z(t) but would not be present in a large- N system. To see if this effect is 
relevant, we estimated the time it would take for an acoustic wave to cross 
the region, using the numbers from [|TJ . The speed of sound in sodium at its 
melting temperature is 2.5 x 10 5 cm/s, and the volume of the region occupied 
by one atom is 278 a 3 ,, so from the fact that there are 500 atoms one finds that 
the time required for an acoustic wave to cross the region is 783 St, or about 
800 timesteps. (The speed of sound in sodium at our lower temperatures 
varies from that at the melting point by roughly 5%, so this result is valid 
to the same accuracy.) In the Figure below, the MD result for Z(t) begins 
to show small oscillatory revivals at about this time; we conclude that this 
is a finite size effect, but it does not affect the data before that time. 

In Fig. H, Eq. fll5|) is plotted on top of the MD data for Z(t) of sodium at 
the two temperatures. Although both temperatures compare exceptionally 
well to the harmonic theory, the match is visibly poorer for the lower temper- 
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ature. However, repeated MD runs at the lower temperature revealed that 
overall variations in Z(t) amount to 10~ 2 on average, which is the same order 
as the differences between theory and MD in this Figure. By 500 timesteps 
the theory is slightly out of phase with the MD data, and this small difference 
persists out to more timesteps. 



4 Conclusions 

These results show that the motion of a liquid in a single potential valley 
is harmonic to an extremely high approximation; the harmonic prediction 
for the function Z(t) matches the calculation from MD very closely. Any 
contributions due to anharmonicity (which are certainly present) are at most 
of the same order as the accuracy of the MD calculations. 

Some form of harmonic approximation, such as the one used here, has 
been taken up by many workers attempting to understand the dynamics of 
liquids, and it is helpful to compare their models with our approach. One of 
the most popular is the theory of Instantaneous Normal Modes (INM) , intro- 
duced by Rahman, Mandell, and McTague and LaViolette and Stillinger 
@ and developed extensively by Stratt (for example, [§]). Stratt expands 
the many-body potential in the neighborhood of an arbitrary point to sec- 
ond order in displacements from that point, and he expresses the potential 
as a quadratic sum of normal modes, in which the frequencies may be either 
real or imaginary. He then replaces the frequencies by their thermal aver- 
ages over the potential surface, resulting in a temperature-dependent density 
of frequencies. From this point he calculates the system's motion and con- 
siders various time correlation functions, including Zit). He observes that 
his results are accurate to order t 4 for short times, but his predictions also 
diverge from MD results very rapidly, in a time shorter than half of one vi- 
brational period. The agreement with MD at long times can be improved 
by omitting the imaginary frequencies from the calculation of Z(t), but of 
course this makes the short time behavior inexact. (The work of Vallauri 
and Bermejo || follows Stratt 's procedure.) Efforts to improve the long time 
behavior of the correlation functions calculated using INM have been made 
by Madan, Keyes, and Seeley |H|, who have attempted to extract from the 
imaginary part of the INM spectrum a damping factor for Z(t) of the general 



type suggested by Zwanzig |ITfl. Also taking their cue from Zwanzig, Cao 



and Voth |T2| have followed a slightly different path, replacing the actual 
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potential by a set of temperature-dependent effective normal modes which, 
as they emphasize, bears little resemblance to the mechanical normal modes 
of a single many-particle valley. In fact, they state quite explicitly that a 
theory based on purely mechanical normal modes will have little success in 
accounting for equilibrium or dynamical properties of liquids. 

An obvious difference between our theory and INM is the nature of our 
approximation. In INM, one approximates the potential quadratically at an 
arbitrary point, with the result that the motion so predicted is accurate only 
for very small times; in our theory, we expand the potential only at very 
special points where we know the predicted motion will be valid for very 
long times. Both theories then face the problem of extending their validity 
beyond the initial approximation, of course, and we will briefly mention our 
extension in the final paragraph below, but there is one particular reason 
why we strongly prefer the approach taken here: The other models all re- 
place the true potential by a temperature-dependent potential determined 
by one or another thermal averaging process. A temperature-dependent po- 
tential does not provide a true Hamiltonian, and therefore it cannot be used 
to calculate the quantum or classical motion, i.e., it cannot be used in the 
Schrodinger equation or Newton's law. (On the dynamical level, tempera- 
ture is not even a well-defined concept.) Further, the Hamiltonian resulting 
from a temperature-dependent potential cannot be used to do statistical me- 
chanics, except through uncontrolled self-consistent procedures. We prefer 
to build our theory in terms of the actual potential, hence in terms of its 
true Hamiltonian, and to find at least approximate solutions for the Hamil- 
tonian motion, so we can apply the standard procedures of equilibrium and 
nonequilibrium statistical mechanics. 

Further, we would argue that Cao and Voth's skepticism regarding purely 
dynamical approaches is unfounded, given the results here. It is difficult to 
compare our Z(t) results to those of others, because their MD-simulated 
states are not always characterized as diffusing or nondiffusing. We are fairly 
confident that Vallauri and Bermejo's Fig. 2b is a comparable state (glassy 
Cs at 20 K), and we believe our fit to MD is slightly better. Madan, Keyes, 
and Seeley's Fig. 3b is an ambiguous case (it is likely that a glass transition 
has occurred), but there also we are confident that our match with MD is 
better. Hence we would claim that this method shows as much promise as 
the others currently available, and with the physical potential as opposed to 
a thermal average potential. 

It is also instructive to compare the results of this paper with a model for 
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Z(t) previously proposed by Wallace |JJ in which a single particle oscillates in 
a three-dimensional harmonic valley, and at each turning point it may with 
probability fi "transit" to an adjacent valley. To apply that model to a non- 
diffusing case, we set /i = (indicating no transits), yielding Z{t) = cos(ut). 
Clearly this would not fit the MD data for any u, and it is easy to see why: 
Wallace included only one frequency in his earlier model, whereas our Eq. 
(P~5| ) contains contributions from many frequencies, all of which are necessary 



to raise the first minimum in Z(t) above —1 and then damp Z(t) out by de- 
phasing. This suggests an alternate path to understanding diffusing states: 
Begin with a mean atom trajectory model that by construction reproduces 
the correct result for Z(t) in the nondiff using regime (Eq. (PD), and then 
incorporate Wallace's notion of transits into this model. Our work in this di- 
rection, with comparison to MD data for higher-temperature diffusing states 
of liquid sodium, will be described in a subsequent paper [T3[ . 
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